%% PuReMD manual

\documentclass{article}

\usepackage{hyperref}
\usepackage{minted}
\usepackage{listings}

\lstset{
  basicstyle=\ttfamily,
  mathescape
}


\title{Manual for PuReMD: \\
  {\bf Pu}ReMD {\bf Re}active {\bf M}olecular {\bf D}ynamics Program}

\author{
  H. Metin Aktulga \\
  \texttt{hma@cse.msu.edu} \\
  \and
  Kurt A. O'Hearn \\
  \texttt{ohearnku@msu.edu}
}

\begin{document}

\maketitle

This manual is for the PuReMD software which has
come to existence as a result of our ReaxFF realization efforts. 
Our initial efforts have led to the SerialReax program, which is a 
sequential implementation for ReaxFF. SerialReax has helped us in verifying 
the accuracy of our implementation in C against the original ReaxFF code 
which was developed in Fortran, such a task would be cumbersome in a parallel 
context. It has also served as a test bed for quickly implementing new 
algorithms and numerical techniques and benchmarking their effectiveness 
before we incorporated them into the parallel version.

PuReMD (PuReMD Reactive Molecular Dynamics) program is based on the 
sequential implemention, SerialReax, and highly efficient and scalable
parallelization techniques. It inherits the excellent 
performance and small memory foot-print features of SerialReax and 
extends these desirable capabilities to systems of sizes that are of
interest to computational scientists. 

For reasons described above, setting up a simulation and running it using 
PuReMD or SerialReax is quite similar to each other. Therefore in this 
manual, we take PuReMD as our basis and describe it first. In a following 
section, we describe the extras that come with SerialReax which we hope 
to incorporate into PuReMD in the near future.


\section{Input Files}
\label{sec:puremd_inp}

PuReMD expects 3 input files: a geometry file describing the system to be 
simulated, a force field file containing ReaxFF parameters and a control 
file to manage simulation variables.


\subsection{Geometry File}
\label{sec:puremd_geo}

Geometry file tells about the types and initial positions of the atoms 
in the system. PuReMD supports geometry files in two formats: 
the PDB format and our custom input format. It is also possible to 
restart from an earlier simulation check-point using a restart file
(which can be either in ASCII or binary format as explained below). 

\subsubsection{PDB format}
\label{sec:puremd_pdb}

For detailed and up-to-date information on the PDB format, please visit 
\href{http://www.wwpdb.org/docs.html}{here}. Input files of various other formats 
can easily be converted to the pdb format using the freely available 
\href{http://openbabel.sourceforge.net/wiki/Main_Page}{OpenBabel software}).

In the geometry file, each atom is assigned a unique serial id to be
able to identify atoms easily during the simulation. PDB format limits the 
number of digits in the atom serial field to only 5 digits, therefore the 
maximum number of atoms that can be input using the PDB format is only 
$100000$, way behind what is generally used for parallel MD simulations.

\subsubsection{custom format}
\label{sec:puremd_custom}

PuReMD features a very simple custom geometry format to alleviate the 
maximum number of atoms limitation of the PDB format and to ease the task of
preparing a geometry file. The general layout of our custom GEO format is as
follows: The first line describes the simulation box and the second line 
gives the total number of atoms in the system. These initial two lines need 
to be followed by a single line for each atom describing it in detail. 
Here is what a custom geo file looks like:
\begin{lstlisting}
BOXGEO x_len y_len z_len alpha beta gamma
N
1 ele1  name1  x1 y1 z1
2 ele2  name2  x2 y2 z2
$\vdots$
N eleN  nameN  xN yN zN 
\end{lstlisting}

First three floating point numbers on the first line give the length of 
the simulation box in x, y, z dimensions, the remaining ones are for the 
angles between each box dimension. Currently, PuReMD works only with an 
orthogonal box meaning all angles need to be $90.0$ degrees. 

There is no limit by the format on the number of atoms that can be input, 
so N shown on the second line can be as large as allowed by the memory
limitations. Each atom line starting from line 3 and until line $N+2$ consists 
of 6 fields: 
\begin{itemize}
  \item an integer denoting the atom's serial id
  \item a string for the chemical symbol of the element 
    (2 characters max, case insensitive)
  \item a string for the atom name (7 characters max)
  \item 3 floating point numbers describing the position in cartesian 
    coordinates
\end{itemize}


\subsection{Force Field File}
\label{sec:puremd_ffield}

Force field file contains the ReaxFF parameters to be used during 
the simulation. Adri van Duin is the main developer and distributor 
for Reax force fields, you can see his contact info
\href{http://www.mne.psu.edu/vanduin}{here}.


\subsection{Control File} 
\label{sec:puremd_control}

Parameters in the control file allow the user to tune various simulation 
options. Parameter names are case-sensitive but their order is not important 
(except that {\tt ensemble\_type} needs to precede both {\tt p\_mass} and 
{\tt pressure}). The table below lists all parameters names which the
software recognizes, and each parameter is described in further detail
below. If a parameter is missing from the control file, its default 
value (as given in each parameter's description below) will be assumed.
Each parameter must be specified in a single line, first token should be
the parameter and the second token should be an appropriate value. 
Comments regarding a parameter can be included after the value field 
on the same line.

\begin{center}
\begin{tabular}{|c|c|} \hline
  \hyperref[sec:simulation_name]{simulation\_name} & \hyperref[sec:ensemble_type]{ensemble\_type} \\ \hline
  \hyperref[sec:nsteps]{nsteps} & \hyperref[sec:dt]{dt} \\ \hline
  \hyperref[sec:reposition_atoms]{reposition\_atoms} & \hyperref[sec:restrict_bonds]{restrict\_bonds} \\ \hline
  \hyperref[sec:tabulate_long_range]{tabulate\_long\_range} & \hyperref[sec:energy_update_freq]{energy\_update\_freq} \\ \hline
  \hyperref[sec:remove_CoM_vel]{remove\_CoM\_vel} & \hyperref[sec:nbrhood_cutoff]{nbrhood\_cutoff} \\ \hline
  \hyperref[sec:bond_graph_cutoff]{bond\_graph\_cutoff} & \hyperref[sec:thb_cutoff]{thb\_cutoff} \\ \hline
  \hyperref[sec:hbond_cutoff]{hbond\_cutoff} & \hyperref[sec:charge_method]{charge\_method} \\ \hline
  \hyperref[sec:cm_q_net]{cm\_q\_net} & \hyperref[sec:cm_solver_type]{cm\_solver\_type} \\ \hline
  \hyperref[sec:cm_solver_max_iters]{cm\_solver\_max\_iters} & \hyperref[sec:cm_solver_restart]{cm\_solver\_restart} \\ \hline
  \hyperref[sec:cm_solver_q_err]{cm\_solver\_q\_err} & \hyperref[sec:cm_domain_sparsity]{cm\_domain\_sparsity} \\ \hline
  \hyperref[sec:cm_solver_pre_comp_type]{cm\_solver\_pre\_comp\_type} & \hyperref[sec:cm_solver_pre_comp_refactor]{cm\_solver\_pre\_comp\_refactor} \\ \hline
  \hyperref[sec:cm_solver_pre_comp_droptol]{cm\_solver\_pre\_comp\_droptol} & \hyperref[sec:cm_solver_pre_comp_sweeps]{cm\_solver\_pre\_comp\_sweeps} \\ \hline
  \hyperref[sec:cm_solver_pre_app_type]{cm\_solver\_pre\_app\_type} & \hyperref[sec:cm_solver_pre_app_jacobi_iters]{cm\_solver\_pre\_app\_jacobi\_iters} \\ \hline
  \hyperref[sec:temp_init]{temp\_init} & \hyperref[sec:temp_final]{temp\_final} \\ \hline
  \hyperref[sec:t_mass]{t\_mass} & \hyperref[sec:t_mode]{t\_mode} \\ \hline
  \hyperref[sec:t_rate]{t\_rate} & \hyperref[sec:t_freq]{t\_freq} \\ \hline
  \hyperref[sec:pressure]{pressure} & \hyperref[sec:p_mass]{p\_mass} \\ \hline
  \hyperref[sec:compress]{compress} & \hyperref[sec:press_mode]{press\_mode} \\ \hline
  \hyperref[sec:geo_format]{geo\_format} & \hyperref[sec:write_freq]{write\_freq} \\ \hline
  \hyperref[sec:traj_compress]{traj\_compress} & \hyperref[sec:traj_format]{traj\_format} \\ \hline
  \hyperref[sec:traj_title]{traj\_title} & \hyperref[sec:atom_info]{atom\_info} \\ \hline
  \hyperref[sec:atom_forces]{atom\_forces} & \hyperref[sec:atom_velocities]{atom\_velocities} \\ \hline
  \hyperref[sec:bond_info]{bond\_info} & \hyperref[sec:angle_info]{angle\_info} \\ \hline
  \hyperref[sec:test_forces]{test\_forces} & \hyperref[sec:molec_anal]{molec\_anal} \\ \hline
  \hyperref[sec:freq_molec_anal]{freq\_molec\_anal} & \hyperref[sec:dipole_anal]{dipole\_anal} \\ \hline
  \hyperref[sec:freq_dipole_anal]{freq\_dipole\_anal} & \hyperref[sec:diffusion_coef]{diffusion\_coef} \\ \hline
  \hyperref[sec:freq_diffusion_coef]{freq\_diffusion\_coef} & \hyperref[sec:restrict_type]{restrict\_type} \\ \hline
  \hyperref[sec:restart_format]{restart\_format} & \hyperref[sec:restart_freq]{restart\_freq} \\ \hline
\end{tabular}
\end{center}

\subsubsection{simulation\_name}
\label{sec:simulation_name}

\begin{verbatim}
  simulation_name    test_puremd
\end{verbatim}
Output files produced by PuReMD will be in 
{\tt simulation\_name.some\_extension} format. Output files will be 
discussed in more detail in Section~\ref{sec:puremd_output}.

Default: {\tt default.sim}

\subsubsection{ensemble\_type}
\label{sec:ensemble_type}

\begin{verbatim}
  ensemble_type    1
\end{verbatim}
{\tt ensemble\_type} denotes the type of the ensemble to be produced by 
PuReMD. Supported ensembles are as follows:
\begin{itemize}
  \item 0: NVE
  \item 1: bNVT: NVT with Berendsen thermostat
  \item 2: nhNVT: NVT with Nose-Hoover thermostat (under testing)
  \item 3: sNPT: semiisotropic NPT with Berendsen's coupling
  \item 4: iNPT: isotropic NPT with Berendsen's coupling
  \item 5: NPT: anisotropic NPT with Parrinello-Rehman coupling 
    (under development)
\end{itemize}

Default: 0 (NVE)

\subsubsection{nsteps}
\label{sec:nsteps}

\begin{verbatim}
  nsteps     1000
\end{verbatim}
{\tt nsteps} controls the total number of steps for the simulation.

Default: 0

\subsubsection{dt}
\label{sec:dt}

\begin{verbatim}
  dt         0.25
\end{verbatim}
{\tt dt} controls the length of each time step (in femtoseconds). 

Default: 0.25

\subsubsection{proc\_by\_dim}
\label{sec:proc_by_dim}

\begin{verbatim}
  proc_by_dim     1 1 3
\end{verbatim}
The distributed memory version of PuReMD uses the
domain decomposition technique to distribute the load
among processors. It currently does not have dynamic load balancing.
{\tt proc\_by\_dim} denotes the desired decomposition of the simulation 
box into subdomains (first integer is the number of equal-length 
partitions in x dimension, second integer is for y dimension and 
the last one is for z dimension). Each subdomain is subsequently assigned 
to a single processor. PuReMD constructs a 3D torus based on the 
{\tt proc\_by\_dim} parameter.

Default: 1 1 1 (single processor)

Note: shared memory versions of PuReMD do not accept this parameter.

\subsubsection{geo\_format}
\label{sec:geo_format}

\begin{verbatim}
  geo_format     0
\end{verbatim}
{\tt geo\_format} parameter informs PuReMD about the format of the 
geometry file to be read. Options are:
\begin{itemize}
  \item 0: custom GEO format
  \item 1: PDB format
  \item 2: ASCII restart file
  \item 3: binary restart file
\end{itemize}
PDB and custom formats were already discussed in Section~\ref{sec:puremd_geo}.
Another option is to resume from an older run by setting {\tt geo\_format}
to 2 (for ASCII restarts) or 3 (for binary restarts) and providing the name 
of the restart file as an argument to PuReMD (instead of the GEO file name).
Then PuReMD will read the box geometry, positions and velocities for all 
atoms in the system from the restart file and continue execution from thereon. 

Default: 0 (custom format)

\subsubsection{restart\_format}
\label{sec:restart_format}
\subsubsection{restart\_freq}
\label{sec:restart_freq}

\begin{verbatim}
  restart_format   1
  restart_freq     0
\end{verbatim}
PuReMD can output restart files in an ASCII format (when 
{\tt restart\_format = 0}) or in a binary format (when {\tt restart\_format = 1}).
While ASCII restarts are good for portability, binary restart files are 
much more compact and does not cause any loss of information due to 
truncation of floating point numbers. Binary restart is the default.

There will not be any restart files output unless {\tt restart\_freq} 
parameter is set to a positive integer. A restart file is named as follows: 
{\tt simulation\_name.resS} where {\tt S} denotes the step that the restart 
file is written.

Defaults: 0 (ASCII), 0

\subsubsection{tabulate\_long\_range}
\label{sec:tabulate_long_range}

\begin{verbatim}
  tabulate_long_range    10000
\end{verbatim}
When set to $m$ (must be a positive integer), {\tt tabulate\_long\_range} 
option turns on the tabulation optimization for computing electrostatics and 
van der Waals interactions. The range [0, cutoff] is sampled at $m$ equally 
spaced points; energy and forces due to long range interactions between each 
atom type in the system are computed at each of these sample points and 
are stored in a table. Then for each interval, coefficients of a fitted
cubic spline interpolation function are computed. During the simulation 
while computing the long range interactions between any two atoms, 
the appropriate interpolation function is located and energy and forces 
between the atom pair is approximated by means of cubic spline interpolation.
This method gives significant speed-up compared to computing everything from 
scratch each time and with only 10000 sample points it is able to provide 
results with an accuracy at machine precision level.

Default: 0 (no tabulation)

\subsubsection{energy\_update\_freq}
\label{sec:energy_update_freq}

\begin{verbatim}
  energy_update_freq     10
\end{verbatim}
This option controls the frequency of writes into output files described 
in detail in Section~\ref{sec:puremd_output} (except for the trajectory 
and restart files which are controlled by other parameters explained
separately).

Default: 0 (no energies in performance logs)

\subsubsection{remove\_CoM\_vel}
\label{sec:remove_CoM_vel}

\begin{verbatim}
  remove_CoM_vel     500
\end{verbatim}
Removal of translational and rotational velocities around the center of 
mass needs to be done for NVT and NPT type ensembles to remove the 
nonphysical effects of scaling velocities. In case of NVE, this is  
unnecessary and is not done regardless of the value of {\tt remove\_CoM\_vel}.

Default: 25

\subsubsection{nbrhood\_cutoff}
\label{sec:nbrhood_cutoff}
\subsubsection{thb\_cutoff}
\label{sec:thb_cutoff}
\subsubsection{hbond\_cutoff}
\label{sec:hbond_cutoff}

\begin{verbatim}
  nbrhood_cutoff     5.0     
  thb_cutoff         0.001   
  hbond_cutoff       7.50
\end{verbatim}
These cutoff parameters are crucial for the correctness and efficiency
of PuReMD. Normally, bonded interactions are truncated after 4-5~\AA\ in 
ReaxFF and this is controlled by the {\tt nbrhood\_cutoff} parameter.

{\tt thb\_cutoff} sets the bond strength threshold for valence angle 
interactions. Bonds which are weaker than {\tt thb\_cutoff} will not 
be included in valence angle interactions.

{\tt hbond\_cutoff} controls the distance between the donor and acceptor 
atoms in a hydrogen bond interaction. Its typical value is from 6\AA\ to 
7.5~\AA. If {\tt hbond\_cutoff} is set to 0, hydrogen bond interactions 
will be turned off completely (could be useful for improved
performance in simulations where it is known apriori that there are no 
hydrogen bonding interactions).

Defaults: 4.0, 0.001, 0.0

\subsubsection{reneighbor}
\label{sec:reneighbor}
\subsubsection{vlist\_buffer}
\label{sec:vlist_buffer}

\begin{verbatim}
  reneighbor     1
  vlist_buffer   2.5
\end{verbatim}
PuReMD features delayed neighbor generation by using Verlet lists. 
{\tt reneighbor} controls the reneighboring frequency and {\tt vlist\_buffer} 
controls the buffer space beyond the maximum ReaxFF interaction cutoff. 

Defaults: 1 (reneighbor every step), 2.5

\subsubsection{charge\_method}
\label{sec:charge_method}
\subsubsection{cm\_q\_net}
\label{sec:cm_q_net}
\subsubsection{cm\_solver\_type}
\label{sec:cm_solver_type}
\subsubsection{cm\_solver\_max\_iters}
\label{sec:cm_solver_max_iters}
\subsubsection{cm\_solver\_restart}
\label{sec:cm_solver_restart}
\subsubsection{cm\_solver\_q\_err}
\label{sec:cm_solver_q_err}
\subsubsection{cm\_solver\_sparsity\_enabled}
\label{sec:cm_solver_sparsity_enabled}
\subsubsection{cm\_solver\_sparsity}
\label{sec:cm_solver_sparsity}
\subsubsection{cm\_solver\_pre\_comp\_type}
\label{sec:cm_solver_pre_comp_type}
\subsubsection{cm\_solver\_pre\_comp\_sweeps}
\label{sec:cm_solver_pre_comp_sweeps}
\subsubsection{cm\_solver\_pre\_comp\_refactor}
\label{sec:cm_solver_pre_comp_refactor}
\subsubsection{cm\_solver\_pre\_comp\_droptol}
\label{sec:cm_solver_pre_comp_droptol}
\subsubsection{cm\_solver\_pre\_app\_type}
\label{sec:cm_solver_pre_app_type}
\subsubsection{cm\_solver\_pre\_app\_jacobi\_iters}
\label{sec:cm_solver_pre_app_jacobi_iters}

\begin{verbatim}
  charge_method                     0
  cm_q_net                          0.0
  cm_solver_type                    0
  cm_solver_max_iters               20
  cm_solver_restart                 100
  cm_solver_q_err                   1e-6
  cm_domain_sparsity                1.0
  cm_solver_pre_comp_type           1
  cm_solver_pre_comp_refactor       1000
  cm_solver_pre_comp_droptol        0.0
  cm_solver_pre_comp_sweeps         3
  cm_solver_pre_app_type            0
  cm_solver_pre_app_jacobi_iters    50
\end{verbatim}
PuReMD uses one of several charge methods for dynamically
determining atomic charges. Unpinning these charge methods
are iterative linear solvers with optional preconditioning.

{\tt charge\_method} controls which charge method is used.
The options are charge equilibration
(0), electronegivity equilibration (1), or atom-condensed Kohn-Sham
approximated to second order (2).

{\tt cm\_q\_net} controls the net system charge.

{\tt cm\_solver\_type} controls which linear solver is used. Options
are GMRES with restarts (0), Householder GMRES with restarts (1),
conjugant gradient (2), and steepest descent (3).

{\tt cm\_solver\_max\_iters} controls the maximum number of iterations
the solver is allowed to perform in order to achieve convergence of
the solution to within the reqiured tolerance.

{\tt cm\_solver\_restart} controls the maximum number of inner iterations
that GMRES-based solvers are allowed before performing a restart.

{\tt cm\_solver\_q\_err} sets the the solution tolerance for the solver.

{\tt cm\_solver\_sparsity} sets the sparsification ratio of the charge matrix
used to compute the preconditioner. Specifically, an additional distance-based
cutoff is applied to compute elements of the (sparsified) charge matrix, and
this matrix is in turn used to compute a preconditioner. The value of this
parameter is multipled by the current neighbor cutoff value to obtain the new
cutoff; hence, a value between 0.0 and 1.0, exclusive, is expected.

{\tt cm\_solver\_pre\_comp\_type} sets the type of preconditioner to be
computed. Options are none (0), Jacobi/diagonal inverse (1), incomplete
Cholesky with dual thresholding (2), incomplete LU computed in an iterative
fashion (3), and iterative incomplete LU single thresholding (4).

{\tt cm\_solver\_pre\_comp\_refactor} sets the number of simulation steps
after which to recompute the preconditioner.

{\tt cm\_solver\_pre\_comp\_droptol} sets the dropping tolerance for
computing incompute Cholesky or LU preconditioners.

{\tt cm\_solver\_pre\_comp\_sweeps} sets the number of sweeps (iterations)
to perform when computing incompute LU iteratively.

{\tt cm\_solver\_pre\_app\_type} determines the type of method used
to apply the preconditioner in the case of two-sided preconditioning
(incomplete Cholesky and LU). Specifically, the application of the
approximate triangular factors requires solving triangular linear systems,
and this parameter controls how this is performed. Options are serial solve
forward/backword substitution (0), solve via level scheduling (1), solve
via graph coloring (3), or approximate solve via Jacobi iteration (4).

{\tt cm\_solver\_pre\_app\_jacobi\_iters} controls how many iterations
are performed for each approximate triangular solve via Jacobi iteration.

{\tt cm\_solver\_freq} can be used to compute charges at every 
few steps instead of the default behaviour of performing it at every 
step. Although doing this less frequently would save important 
computational time, it is not recommended. Because this might cause wild 
fluctuations in energies and forces.

Defaults: 

Notes:
\begin{enumerate}
  \item Only the shared memory versions of PuReMD contain all the solvers,
    while only CG with diagonal inverse preconditioning is contained in the
    distributed memory versions.
  \item The full set of preconditioning options is implemented in the shared
    memory non-GPU version of PuReMD. The other versions contain only diagonal
    inverse preconditioning.
\end{enumerate}

\subsubsection{temp\_init}
\label{sec:temp_init}
\subsubsection{temp\_final}
\label{sec:temp_final}
\subsubsection{t\_mass}
\label{sec:t_mass}

\begin{verbatim}
  temp_init    0.0
  temp_final   300.0
  t_mass       0.1666
\end{verbatim}
Temperature coupling parameters ({\tt temp\_final} and {\tt t\_mass}) are 
effective in all types of ensembles except for NVE. Initial temperature 
is controlled via the {\tt temp\_init} parameter including the NVE ensemble.
0~K is the default value for {\tt temp\_init} and 300~K is the default value 
for {\tt temp\_final}. PuReMD features both Berendsen~\cite{ref:berendsen} 
and Nose-Hoover~\cite{ref:klein} type thermostats as was mentioned while 
explaining the {\tt ensemble\_type} parameter.

{\tt t\_mass} is the thermal inertia given in femtoseconds. Suggested
value of {\tt t\_mass} is 500.0 and 0.166 for the Berendsen 
thermostat and the Nose-Hoover thermostats, respectively.

Defaults: 0.0, 300.0, 0.16666

Note: Nose-Hoover thermostat is still under testing

\subsubsection{pressure}
\label{sec:pressure}
\subsubsection{p\_mass}
\label{sec:p_mass}

\begin{verbatim}
  pressure      0.000101 0.000101 0.000101
  p_mass        5000.0   5000.0   5000.0
\end{verbatim}
Pressure coupling parameters are needed only when working with NPT-type 
ensembles. Currently iNPT (isotropic NPT) and sNPT (semi-isotropic NPT) 
are the available pressure coupling ensembles in PuReMD. Berendsen 
thermostats and barostats are used in both cases~\cite{ref:berendsen}. 
{\tt pressure} is the desired pressure of the system in GPa and {\tt p\_mass}
is the virial inertia in femtoseconds. Suggested (and the default) value of
{\tt p\_mass} is 5000.0 together with a {\tt t\_mass} of 500.0 as NPT methods
use Berendsen-type thermostats only.

For the iNPT ensemble, {\tt pressure} parameter expects a single
floating number (in case there are more, they will simply be ignored) 
to control pressure. For the sNPT ensemble, {\tt pressure} parameter 
expects 3 floating point numbers to control pressure on each dimension.
Same things apply for {\tt p\_mass} as well.

Defaults: 0.000101325, 5000.0

\subsubsection{write\_freq}
\label{sec:write_freq}
\subsubsection{traj\_format}
\label{sec:traj_format}

\begin{verbatim}
  write_freq     100
  traj_format      1
\end{verbatim}
Trajectory of the simulation will be output to the trajectory file 
(which will automatically be named as {\tt simulation\_name.trj}) at 
every {\tt write\_freq} steps. For making analysis easier, the trajectory 
file is written as an ASCII file.

The distributed memory version of PuReMD
can output trajectories either using simple MPI send/receives 
(option 0 which is the default) or using MPI I/O calls (option 1) which 
are part of the MPI-2 standard. The latter option is supposed to be more 
efficient (not verified by tests though) but may not be available in some 
MPI implementations. {\tt traj\_method} option is not applicable to 
SerialReax simulations.

Defaults: 0 (no trajectory file written), 0 (simple I/O)


\subsubsection{traj\_title}
\label{sec:traj_title}
\subsubsection{atom\_info}
\label{sec:atom_info}
\subsubsection{atom\_forces}
\label{sec:atom_forces}
\subsubsection{atom\_velocities}
\label{sec:atom_velocities}
\subsubsection{bond\_info}
\label{sec:bond_info}
\subsubsection{angle\_info}
\label{sec:angle_info}

\begin{verbatim}
  traj_title          TEST
  atom_info           1
  atom_forces         1
  atom_velocities     1
  bond_info           0
  angle_info          0
\end{verbatim}
Currently PuReMD only outputs trajectories in its custom trajectory 
format. This custom format starts with a trajectory header detailing 
the trajectory title and the values of control parameters used for 
the simulation. A brief description of atoms follow the trajectory 
header with atom serial ids and what element each atom is.

Then at each {\tt write\_freq} steps (including step 0), a trajectory 
frame is appended to the trajectory file.  The frame header which gives 
information about various potential energies, temperature, pressure and 
box geometry is standard. However, the latter parts of the frame can be 
customized using {\tt atom\_info}, {\tt atom\_forces}, {\tt atom\_velocities}, 
{\tt bond\_info} and {\tt angle\_info} parameters which are already
self-explanatory. The ordering is atoms section, bonds section and angles 
section assuming that they are all present.

One nice property of the custom trajectory format is that each part of 
the trajectory is prepended by a number that can be used to skip that part.
For example, the trajectory header is prepended by an integer giving the 
number of characters to skip the control parameters section. The initial 
atom descriptions is prepended by the number of characters to skip the 
initial descriptions part and another one that tells the number of atom 
description lines. Similar numbers are found at the start of each section 
within a trajectory frame as well, making it easy to skip parts which are 
not of interest to a particular trajectory analysis procedure. So the 
general layout of our custom trajectory format is as follows (assuming 
all trajectory options are turned on):
\begin{lstlisting}
CHARS_TO_SKIP_SECTION
trajectory header
CHARS_TO_SKIP_ATOM_DESCS NUM_LINES
atom descriptions
CHARS_TO_SKIP_FRAME_HEADER
frame1 header
CHARS_TO_SKIP_ATOM_LINES NUM_ATOM_LINES
frame1 atom info
CHARS_TO_SKIP_BOND_LINES NUM_BOND_LINES
frame1 bond info
CHARS_TO_SKIP_ANGLE_LINES NUM_ANGLE_LINES
frame1 angle info
$\vdots$
CHARS_TO_SKIP_FRAME_HEADER
frameN header
CHARS_TO_SKIP_ATOM_LINES NUM_ATOM_LINES
frameN atom info
CHARS_TO_SKIP_BOND_LINES NUM_BOND_LINES
frameN bond info
CHARS_TO_SKIP_ANGLE_LINES NUM_ANGLE_LINES
frameN angle info
\end{lstlisting}

Defaults: default\_title, 0, 0, 0, 0, 0 (all off)


\section{SerialReax Extras}
\label{sec:serialreax_extras}

In this section, we explain the parameters found in SerialReax but not in
PuReMD. Our work towards adding the same functionalities into PuReMD is 
underway.

%In addition to the PCG solver, SerialReax features a preconditioned GMRES 
%(PGMRES) solver and an incomplete LU factorization (ILU) based  
%preconditioning scheme. An ILU factorization essentially does the 
%same thing as an LU factorization but small terms in the matrix are dropped
%to expedite the factorization and to prevent a huge number of fill-ins in the
%factor matrices. Following are the extra control parameters found in 
%SerialReax regarding the QEq solver:
%\begin{verbatim}
%  ilu_refactor        100
%  ilu_droptol         0.01
%\end{verbatim}
%{\tt ilu\_droptol} sets the threshold for dropping small terms in the 
%resulting ILU factors. Suggested (and the default) value for 
%{\tt ilu\_droptol} is $10^{-2}$. Despite the drop rules, ILU factorization 
%is still a costly operation. So a user can choose to perform it at 
%every {\tt ilu\_refactor} steps. The fact that atoms move very slowly in an 
%MD simulation allows the use of same ILU factors as preconditioners in the 
%subsequent steps with little performance loss. For liquids, this frequency 
%can be on the order of 100-200 steps, for solids it can go up to thousands 
%of steps depending on how fast atoms are moving. The default for 
%{\tt ilu\_refactor} is 100.

\begin{verbatim}
  t_mode        0
  t_rate     -100.0
  t_freq        2.0
\end{verbatim}
These options  are specifically for being able to change the temperature of 
the system during a simulation. {\tt t\_mode} of 1 gives a step-wise 
control over temperature, \emph{i.e.} the system maintains its temperature
for a simulation time of {\tt t\_freq} picoseconds. After that, the 
target temperature is increased/decreased by {\tt t\_rate}~K and the 
thermostat lets the system converge to its new target temperature.

On the other hand, {\tt t\_mode} of 2 increases/decreases the target 
temperature at each time-step by an amount which corresponds to 
{\tt t\_rate / t\_freq}~K/ps. The overall effect of such a regime is a 
constant slope (instead of the step pattern with a {\tt t\_mode} of 1) 
in the target temperature--simulation time graph.

\begin{verbatim}
  molec_anal          1
  freq_molec_anal     1
  bond_graph_cutoff   0.3
  ignore              2  0 3
\end{verbatim}
Since ReaxFF is a reactive force field, during the simulation molecules 
present in the system will change. These changes can  be traced by turning 
the {\tt molec\_anal} option on by setting it to a non-zero integer. 
Molecules are determined based on the {\tt bond\_graph\_cutoff} parameter: 
bond orders less than this threshold are not counted as a physical
bond and do not contribute to molecular structures, all others do. 
{\tt ignore} allows one to ignore the bondings of specific atom types. 
The first number after {\tt ignore} denotes how many atom types will be 
listed and the following numbers (which correspond to the order of 
elements in the force field file) denote the atom types to be ignored 
in molecular analysis.


%\begin{verbatim}
%  dipole_anal         1    ! 1: calc electric dipole moment
%  freq_dipole_anal    1    ! electric dipole calc freq
%\end{verbatim}
%Currently electric dipole moment analysis is available for water molecules 
%only but it can be extended to other molecule types upon request.


%\begin{verbatim}
%  diffusion_coef      1    ! 1: calc diffusion coef
%  freq_diffusion_coef 1    ! diffusion coef calc freq
%  restrict_type       2    ! -1: all, >=0: only this type
%\end{verbatim}
%Our program allows you to compute the diffusion coefficient of 
%the system, too. It can be restricted to certain types of atoms 
%if desired.

%\begin{verbatim}
%  reposition_atoms        0    ! 1: CoM-centered, 2: CoM-origin
%\end{verbatim}
%This option lets you position the atoms based on your choice. 
%Option $0$ will simply fit the atoms into the periodic box. 
%Option $1$ will translate the atoms such that the center of mass 
%will be at the center of the box, 
%whereas option $2$ positions the center of mass to the origin. 
%Options $1$ and $2$ may need further testing, so it is safer to 
%use option $0$ for now.

%\begin{verbatim}
%  restrict_bonds          0    ! turns reactions on and off
%\end{verbatim}
%When set to $m$ (must be a positive integer), this option enforces 
%bonds given in CONECT lines of geometry file for the first $m$ 
%steps of the simulation. 
%This is done by including only the atoms given on the CONECT lines 
%among the near neighbors of an atom. 
%Turning this option would probably produce nonphysical trajectories 
%but may be useful for energy minimization purposes.

% \begin{verbatim}
%  periodic_boundaries  1       ! 1: periodic boundaries on 
%  periodic_images      3 3 3   ! no of images in each direction
% \end{verbatim}
% Above parameters are concerned with the boundary conditions of the
% system. 
% Currently only periodic boundaries are supported but non-periodic boundaries 
% effect can be accomplished by making the simulation box big enough. 
% Note that adding an empty space of \emph{r\_cut} \AA (discussed below) in 
% \emph{x,y,z} dimensions will be enough for a completely non-periodic 
% simulation box. 
% Or if desired, any combination of dimensions might be made non-periodic
% by adding this empty space to them and letting others end where the system 
% ends.
% Since currently we are implementing shielded electrostatic interactions, 
% \emph{periodic\_images} is not effective yet. It will be required when 
% electrostatic interactions across periodic images are implemented.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Compilation and Execution}
\label{sec:puremd_execute}

PuReMD is distributed in the \mintinline{bash}{tar.gz} compression format which can 
be extracted under a Unix system with the following command:
\mint{bash}{tar -xvf PuReMD.tar.gz}

This results in the creation of a new directory, named \mintinline{bash}{PuReMD}, which will
appear in the working directory. The base directory contains a configure script
along with the necessary makefiles to compile a particular version of the
software; these steps are elaborated further below.  (For developers only: see
the documentation on the Gitlab server for how to generate these files using
the GNU Autotools) In addition, there is also a directory with several sample
systems and force field files (\mintinline{bash}{data/benchmarks}), and a
directory with relevant documentation including this user manual
({\mintinline{bash}{doc}).

To build PuReMD, you must run the following sequence of commands:
\begin{minted}{bash}
  ./configure
  make
  make install # optional
\end{minted}

Upon successful compilation, a binary executable file will be located
in \mintinline{bash}{*/bin}, where \mintinline{bash}{*} is one of the
following directories:
\begin{minted}{bash}
  sPuReMD    # serial/OpenMP shared memory code
  PuReMD-GPU # GPU shared memory code
  PuReMD     # MPI distributed memory code
  PG-PuReMD  # MPI+GPU shared memory code
\end{minted}

There are several different versions of PuReMD which can be compiled,
and these versions can be enabled or disabled via the following options
to the configure script:
\begin{minted}{bash}
  --enable-serial=yes  # build serial shared memory version
  --enable-openmp=yes  # build OpenMP shared memory version
  --enable-gpu=yes     # build GPU shared memory version
  --enable-mpi=yes     # build MPI distributed memory version
  --enable-mpi-gpu=yes # build MPI+GPU distributed memory version
\end{minted}
Furthermore, run \mintinline{bash}{./configure --help} to see the full list
of variables and options which can be set.

PuReMD requires 3 input files as mentioned in section~\ref{sec:puremd_inp}. 
For example, the command to run \mintinline{bash}{puremd} with OpenMPI is as follows:
\begin{minted}{bash}
  mpirun -np num_procs -machinefile m.txt PuReMD/bin/puremd
    path/to/geo path/to/ffield path/to/control
\end{minted}

SerialReax comes in a similar distribution format and Makefile,
so instructions for compiling and running PuReMD is applicable for 
SerialReax as well.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Output}
\label{sec:puremd_output}

PuReMD writes its output files into the directory where it is run. 
There are a number of output files all of which have the 
{\tt simulation\_name} as the first part of their names followed 
by its unique extension:

\begin{itemize}
  \item{\textbf{.out}} contains a summary of the simulation progress. 
    Its format is:
    \begin{verbatim}
    Step   Total Energy   Potential   Kinetic   
    T (in K)    Volume(in A^3)       P(in GP)
  \end{verbatim}

\item{\textbf{.pot}} contains detailed information regarding various types of
  energies that comprise the total potential energy:
    \begin{verbatim}
    Step   Bonds   OverCoor+UnderCoor   LonePair   
    Angle+Penalty   3-body Coalition    Hydrogen Bonds  
    Torsion   4-body Conjugation   
    vander Waals   Coulomb   Polarization
  \end{verbatim}

\item{\textbf{.log}} is intended for performance tracking purposes. 
  It displays the total time per step and what parts of code take 
    up how much time to compute.

  \item{\textbf{.prs}} is output only when pressure coupling is on. 
    It displays detailed information regarding the pressure and 
    box dimensions as simulation progresses.

  \item{\textbf{.trj}} is the trajectory file. Atom positions are written 
    into this file at every {\tt write\_freq} steps using the desired format 
    as explained before. Each frame is concatenated one below the other.

    %\item{\textbf{.dpl}} is for dipole moment analysis:
    %\begin{verbatim}
    %Step   Molecule Count   Avg Dipole Moment
    %\end{verbatim}
    %      
    %\item{\textbf{.drft}} is for diffusion coefficient analysis:
    %\begin{verbatim}
    %Step   Type     Count   Avg Squared Displacement
    %\end{verbatim}
\end{itemize}

Apart from these, there might be some text printed to \emph{stderr} 
for debugging purposes. If you encounter some problems with the code
(like a segmentation fault or unexpected termination of the code),
please contact \href{mailto:hma@cse.msu.edu}{here} with the error message 
printed to \emph{stderr} and your input files.

In addition to the output files above, SerialReax can output another
file (with extension \textbf{.mol}) which contains the fragmentation 
analysis output.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{thebibliography}{99}
  \bibitem{ref:klein}
    Glenn J. Martyna, Douglas J. Tobias, and Michael L. Klein. 
    ``Constant pressure molecular dynamics algorithms.'' 
    The Journal of Chemical Physics 101, 4177 (1994).

  \bibitem{ref:berendsen}
    H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and 
    J. R. Haak.
    ``Molecular dynamics with coupling to an external bath.''
    The Journal of Chemical Physics 81, 3684-3690 (1984).

\end{thebibliography}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}
